%% Posterior function

function PS = posterior_fun_MA(xx,ydata,prior,NS,Nx,NC,tr,setmin)

lnL_sbl_xC = 189*ones(NC,1);
xx = xx';

if tr == 1  % --> parameters in R, transform it to right domain
    alphap11 = exp(xx(1))/(1+exp(xx(1))); thch.p11 = alphap11*prior.p11lb + (1-alphap11)*prior.p11ub;
    alphap22 = exp(xx(2))/(1+exp(xx(2))); thch.p22 = alphap22*prior.p22lb + (1-alphap22)*prior.p22ub;
    alphamu1 = exp(xx(3))/(1+exp(xx(3))); thch.mu1 = alphamu1*prior.mu1lb + (1-alphamu1)*prior.mu1ub;
    alphaDmu = exp(xx(4))/(1+exp(xx(4))); thch.Dmu = alphaDmu*prior.Dmulb + (1-alphaDmu)*prior.Dmuub;   
    alphasg1 = exp(xx(5))/(1+exp(xx(5))); thch.sg1 = alphasg1*prior.sg1lb + (1-alphasg1)*prior.sg1ub;    
    alphasg2 = exp(xx(6))/(1+exp(xx(6))); thch.sg2 = alphasg2*prior.sg2lb + (1-alphasg2)*prior.sg2ub;    
    
    thch.mu2 = thch.mu1 + thch.Dmu;
else
    thch.p11 = xx(1);
    thch.p22 = xx(2);
    thch.mu1 = xx(3);
    thch.Dmu = xx(4);
    thch.sg1 = xx(5);    
    thch.sg2 = xx(6);    
    
    thch.mu2 = thch.mu1 + thch.Dmu;
end

thch.Pmat = [thch.p11 ,  1-thch.p11;...
             1-thch.p22,  thch.p22];
         
thch.R = 0;
thch.Q = 1;

thch.A = [0, 0 ; ...
          1, 0];
           
thch.G = [1; 0];

thch.muvec = [thch.mu1, thch.mu2] ;
thch.sgvec = [thch.sg1, thch.sg2]; 

thch.x00 = zeros(Nx,1);
thch.P00 = eye(Nx);

% run Compute_XI_ss
Compute_XIss
thch.XI_00 = XI_ss;

PP = 1;

for ii = 1:NC
    yaux = ydata(:,ii);
    in   = find(isnan(yaux) == 0, 1, 'first');
    fn   = find(isnan(yaux) == 0, 1, 'last');
    y    = yaux(in:fn);
    
    [~, ~, ~,~, ~, lnL_sbl_xC(ii)] = Kim_Filter_MA_sg(thch,y,NS,Nx);    
    lnL_sbl_xC(ii)                 = log(PP) + lnL_sbl_xC(ii) ;    
    
end

PS = sum(lnL_sbl_xC);

if setmin == 1
    PS = - PS;
end

return